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Abstract 

A rigorous mathematical framework for analyzing the chemical master equa- 
tion (CME) with bistability, based on the theory of large deviation, is proposed. 
Using a simple phosphorylation-dephosphorylation cycle with feedback as an exam- 
ple, we show that a nonequilibrium steady-state (NESS) phase transition occurs in 
the system which has all the characteristics of classic equilibrium phase transition: 
Maxwell construction, discontinuous fraction of phosphorylation as a function of the 
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kinase activity, and Lee- Yang's zero for the generating function. The cusp in non- 
hnear bifurcation theory matches the tricritical point of the phase transition. The 
mathematical analysis suggests three distinct time scales, and related mathematical 
descriptions, of (i) molecular signaling, (ii) biochemical network dynamics, and (iii) 
cellular evolution. The (i) and (iii) are stochastic while (ii) is deterministic. 



Microscopic, stochastic molecular fluctuations disappear in the thermodynamic 
limit in which deterministic nonlinear behavior arises. However, in the mesoscopic 
world of cellular biology, complex dynamics with multiple time scales makes the 
meaning of thermodynamic limit only relative of one level of description with respect 
to another. More specifically, we shall show in the present paper that there are three 
biologically significant time scales, with related levels of mathematical description: 
(i) stochastic molecular signaling, (ii) deterministic biochemical network dynamics, 
and finally (iii) stochastic (again!) cellular evolution. In other words, there is 
stochastic behavior beyond the deterministic dynamics; all three levels are contained 
in a mesoscopic living system. Current cellular molecular biology chiefiy focuses on 
(i) while increasingly interested in (ii); however, it is the (iii), we believe, that is 
most relevant to major cellular biological issues such as differentiation, apoptosis, 
and cancer immunoediting. 

Our conclusion is reached through a detailed mathematical analysis of a simple 
cellular signaling module: a phosphorylation-dephosphorylation cycle (PdPC) with 
feedback p. We use the chemical master equation (CME) as the starting model. 
In recent years CME has emerged as one of the physiochemical foundations of cel- 
lular biochemistry [2]. The theory had begun in 1940 and went through a major 
development in 1960s and 70s [3]. In particular, the Brussels school has used this 
theory as a mathematical bases of nonequilibrium steady state (NESS), a term first 
proposed by Klein [1]. It is now widely accepted that both concepts of CME and 
NESS are appropriate for studying isothermal, homeostatic cellular biochemistry 
[5]. The mathematical theory of NESS is an irreversible, but stationary stochastic 
processes, associated with which the concepts of entropy production and stationary 
distribution naturally arise [6]. 

It is generally believed that the deterministic nonlinear dynamics, derived from 
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the CME in the hmit of the reaction system volume ^ oo according to Kurtz's 
theory [7j, is the macroscopic counterpart of the chemical reaction system [2]. While 
this is certainly true, here we refine this notion by studying the large-deviation 
properties ofV^ oo, i.e., the thermodynamic limit. We shall show that in the case 
of a nonlinear dynamical system with multiple dynamic attractors, there is a unique 
macroscopic thermodynamic state; all the other macroscopic attractors are in fact 
metastable, with an infinitesimal stationary probability oc e~^^ and exponential 
small exit rate oc e~"^ > 0). 

The mathematical theory of large deviation (LDT) [8] is the natural device for 
understanding the thermodynamic limit of systems with mutlist ability, i.e., phase 
transition(s) [9]. Our result based on the LDT rigorously establishes the stochas- 
tic dynamics with bi(multi)-modal distribution as the mesoscopic signature of a 
nonlinear dynamics with bi(multi)-stability. 

We have recently re-examined the nonlinear bistability in the context of biochem- 
ical signaling module [10]. In the thermodynamic limit when V tends infinity, there 
is a phase transition associated with the conventional nonlinear dynamic approach 
based on the Law of Mass Action, which is the macroscopic limit, in some sense, of 
the CME [7]. A Maxwell-type construction is an integral part of a complete theory 
of the CME dOj. 

In equilibrium phase transition, Lee- Yang theorem for grand canonical partition 
function is widely considered to be a deep and elegant result [12]. We shall show, 
non-differentiability of a function c(A) (the NESS counterpart of the free energy 
function) is the origin of multi-phase behavior, and it is because a zero of G{X) (the 
NESS counterpart of a partition function) reaches the real axis. Different attempts 
have been made to generalize the Lee- Yang theory to NESS and to bimodalities in 

m- 

Large deviation theory, Maxwell construction and first-order phase 
transition in a NESS. We now consider the same biochemical signaling system in 
[H [To] in terms of a one-dimensional CME. Let pv{n) as its stationary probability 
for Ny, the random variable representing the activated kinase molecule X; V being 
the volume of the system. 

According to the classic result of LDT [8], [9], especially Sec. 4.5.2 in the text by 
Dembo and Zeitouni's text, it concludes that if ^ satisfies the LDT with a good 



3 



"rate function" i.e., pv{n) ~ e^^'^^^\ x > 0, then 

(a) For each A, the "free energy function" c(A) = hniy^oo y^og{e^'^^) exists, and it 
is finite and nondecreasing. Moreover it satisfies 

c(A) = supjAa; — (p{x)}. (1) 

x>0 

(b) If (f){x) is convex, then it is the Fenchel-Legendre transform of c(A), namely, 

(f){x) = c*{x) = sup{Ax — c(A)}. (2) 
xen 

(c) If is not convex, then c*{x) is the affine regularization of i.e. c*(-) < 
4>{-), and for any convex rate function / such that /(•) < 0(-) implies /(■) < c*(-). 

Consequently, we know that when (j){x) is bimodal with two local minima, then 
they are at different heights if and only if the c(A) is differentiable at A = according 
to the well-known Gartner- Ellis theorem [9], and ^^^^ is simply the position of the 
lower minimum. This implies that the Maxwell construction corresponds to the 
function c(A) being non-analytic at A = 0. Further, if the rate function is 
analytic, then c(A) is continuous and 

c(A) = sup {Xx — (3) 

<j>'{x)=X 

If a non-convex (j){x) has two local minima of 0(a;) with equal height, for suf- 
ficiently small A < 0, c(A) = Xx — (j){x) for the x near the left minima satisfying 
(j)'{x) = A; and when A > 0, also c(A) = Aa; — (j){x) for the x near the right minima 
satisfying (f)' (x) = X. Therefore, the left and right derivatives of the function c(A) at 
A = both exist but are equal to the left and right local minima respectively. 

If (f){x) has two local minima Xi and X2 with different heights, one can rewrite 
Ax — = X'x — ipi^x) where iplx) = (f){x) — X*x such that ip{x) has two minima 
with equal heights and A' = A — A*. Hence, the nonanalytic point of c(A) moves 
to A*, and also the left and right derivatives at A = A* both exist and equal to the 
left and right local minima of tplx) respectively. In other words. A* is just the slope 
of the tangent line of (pix) with exactly two tangent points. More generally, if the 
non-convex 0(x) has k tangent lines with more than one tangent points, then the 
function c(A) has k non-analytical points, and vise versa. So it is a "higher level" of 
convexity! [11] 
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Figure 1: Bimodal, non-convex (f){x) in (A) gives rise to non-analyticity in the 
Fenchel-Legendre transform c(A) in (B), the generating function of the Ny in the 
thermodynamic hmit. In equihbiurm statistical mechanics, c(A) is the free energy; 
hence the non-differentiabihty at A = A* indicates first-order phase transition. A 
quadratic function, with curvature cr~^ in (A) gives a quandratic function, with 
curvature cr^, in (B). c*{x) in (C) is the Fenchel-Legendre transform of c(A), also 
known as the affine regularization of 



The above LDT results are summarized in Fig. [H Fig. [T] shows that, as that in 
Lee- Yang's theory [12], c(A) is continuous but non-differential at A = A*. 

Now let us consider another parameter 6 of the system. Let it be a bifurcation 
parameter in the nonlinear dynamics according to the Law of Mass Action [1] . We 
have shown in that the stable and unstable fixed points of the nonlinear dynam- 
ics correspond precisely with the minima and maxima of the (pix), and bistability 
corresponds to double- wells in and bimodality of —(f){x). 

Here consider the function (l/l^) log(e^^^) = cy(A,^). As V tends to infinity, 
the limit c(A, 6) exists, and it is continuous and a non-decreasing function of A. 
Furthermore, there is a line in the (A, 6) plane at which the c is non-differentiable 
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with respect to A. The hne passes (0,6'*) where 6* is the critical value of Maxwell 



□ 




Figure 2: (A) The solid line represents the extrema of the which corresponds 
to the stable and unstable fixed points of the nonlinear differential equation model. 
6' is a bifurcation parameter. When 9 = 6*, the two wells of (pix) have equal height. 
(B) For each value of 6, the double- well yields a non-analytical point X*{6). 
This line crosses the A = when 6 = 6*. (C) as a function of (D) is in 
fact the position of the lower minima of which is the mean concentration of X 
in the system. 

In our theory, the derivative at A = is particularly meaningful: It is the mean 
concentration of molecules in the system (property of the generating function). In 
the thermodynamic limit, the mean and the highest peak position of e"^"^*^^^ are the 
same, the macroscopic value. Thus, we understand that the Maxwell construction 
implies the mean concentration is not continuous. 
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Generalizing Lee- Yang's theory. In equilibrium phase transition, according 
to [12], the non-analyticity in the free energy function c(A) is due to a zero in the 
partition function Gv{^) = (e^^^) approaching the real axis from the complex plane 
of A. Is the non-analyticity in our c(A) also due to the zero of G{\) = limy^oo G'v(A)? 
This is indeed the case. 

Our probability distribution for Ny has a finite support. So the generating 
function is a finite order polynomial of z, (use z = e^). Then consider a region 
of the complex plane of z, which contains a section of the z axis. According to 
Theorem 2 in [12] which is a pure mathematical result, the zero of the generating 
function must be "pinched" onto real z-axis at the non-analytic point of the free 
energy function c(A) when V tends to infinity. Therefore, our theory generalizes the 
Lee- Yang theory to nonequilibrium phase transition. 

Several previous works have generalized the Lee- Yang theory in nonequilibrium 
steady states [I3] through specific examples. It has been suggested that the bimodal 
distribution could imply the Lee- Yang theory, but not vice versa. This is consistent 
with our result. 

Cusp catastrophe and tri-critical point in a PdPC with feedback. We 

consider the simple PdPC with positive feedback which exhibits nonlinear bistability 

E + K*^E* + K*, K + 2E*^K*, E* + P ^ E + P, (4) 

in which K and K* are inactive and active forms of a kinase, P is a phosphatase. 
E* is the phosphorylated E, a signaling molecule. Usually E* is functionally ac- 
tive, i.e., "turned-on". Following the previous treatment [U HO], we assume the 
reversible binding K + 2E* ^ K* is rapid. Hence, the dynamics of the fraction of 
phosphorylated E, x, satisfies 

dx 

— = Ox'^ [(1 — x) — ex] + [/i(l — x) — x] = r(x; 6, e), (5) 

in which the three parameters 9 represents the ratio of the activity of the kinase to 
that of the phosphatase; e represents the ADP to ATP concentration ratio, and 
represents the strength of phosphorolysis. — /cBTln(/ie) = AG represents the ATP 
hydrolysis energy. In a living cell, both u and e are small; hence 7 = — ^ 1. 

For large system's volume V , the CME gives the stationary probability p"'^'^^(x) oc 



7 



e ^<^(^)^ where the LDT rate function [TU] 

[1 - x){ex^ + 12) 



0(a;) = ln(l-a;)-xln 

x{Oex^ + Ij 
One can easily check that 



/i 




+2^ / — arctan 4 / —x — arctan v Bex. 



(6) 



_ , (l-x)(gx^ + /i) 
rfx a;(^ex2 + 1) ' ^ ^ 

and the extrema match exactly with the roots of r(x; ^, e) = 0. 

The Eq. ^ exhibits saddle-node bifurcations and cusp catastrophe. One obtains 
the parameter region for the bistability from simultaneously solving r(z) =0 and 

'^''(z) _ g. 
dz 

Oz^ [1 - (1 + e)z] + [/i - (1 + fi)z] = 0, e [2z - 3(1 + e)z^] - (1 + /x) = 0. (8) 

The two equations give the boundary of the region of bistability in {6, e) space (in 
terms of z as a parametric curve): 

,^2(l+^_3^ 2,-i, + l)z _ 

z - 2(/i + 1)^2 V ; 

Fig. [3JA. shows the steady states of Eq. (JSj), function of 6 with various e. 

We see for the range of e < 1.33 the system has three fixed points, i.e., bistability. 
After introducing the Maxwell construction for each and every curve x^^{6), we 
obtain a set of monotonic x'^'^{6). This corresponds to the "PV-isotherm" in the van 
der Waals theory of phase transition. 

According to the cusp catastrope theory [T3] , there is a region in the e — 6 plane 
with three fixed points. The boundary of the region is where the system has exactly 
two fixed points, i.e., where bifurcation occurs: ^i(e) and 6*2 (e). One of the most 
important features of this region is that it has a cusp, at 9cusp = ^"^3^^ , Gcusp = 
when Zcusp = as shown in Fig. [3B- 

For a given e, the critical 6* at which the Maxwell construction is performed 
satisfies Oi < 9* < 62. Thus, the critical line 9*{e) abruptly terminates at the cusp. 
In equilibrirum phase transition, the cusp is also known as tri-critical point |15j . 
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We also note that bistability implies that the x'*** as function of the 6, or e, is not 
monotonic (it is S-shaped). However, after the Maxwell construction the resulting 
x'^^{6) is monotonic in "true" thermodynamic limit! It is precisely the same situation 
as the PV isotherm in gas-liquid phase transition. The word "true" means one has 
to wait sufficiently long to allow the jumps back and forth between attractors. The 
biological significance of monotonicity remained to be elucidated. 




Figure 3: (A) Steady state a;** as functions of 6 according to Eq. ([5]), together with 
Maxwell constructions under the shaded region. The parameter used fi = 0.05, with 
e = 1.3, 1,0.5,0.005 from top to bottom. (B) The solid lines represent the saddle- 
node critical points, i.e., the filled circles in (A). They meet at a cusp. The dashed 
line represents the critical value at which the Maxwell construction is performed. 
The dashed line terminates at the cusp. 

Discussion. The present paper shows that many classical concepts from equi- 
librium phase transition can be applied to bifurcation problem in nonlinear chemical 
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djTiamics that has a mesoscopic stochastic underhne in terms of the CME. From the 
CME point of view, the LDT treatment we present is a small but significant step 
beyond the Kurtz theory towards the macroscopic nonlinear dynamics. Analyzing 
the CME is a much more challenging problem than analyzing the partition function 
since the former offers a dynamic theory. 

The celebrated Maxwell construction is a natural consequence of the general 
theory we propose, and the well-known Lee- Yang theorem is in fact a special case 
of it. More importantly, the general theory is applicable to driven systems with 
nonequilibrium steady state. 

On the mathematical side, the general theory provides a framework to study 
nonlinear bifurcations in terms of mathematical non-analyticity of a certain function, 
a vision long being hold by some investigators p^. The large deviation function 
can be in fact considered as some type of stochastic landscape (potential, Lyapunov 
function in a not rigorous sense) for systems without gradient, nor detailed balance 

While the CME as a fundamental theory of studying cellular biochemistry re- 
mains to be validated experimentally, it is certainly an acceptable mathematical 
model for studying mesoscopic complexity and emergent organization, as called by 
Laughlin et al. [TH]. Chemical reactions are marvellous systems for understanding 
complexity. The present work shows that while Kurtz's theorem is correct, the real 
limit of V tends infinity is not the solution to the law of mass action, but rather 
requires a LDT treatment. 

The existence of "nice" (j){x) in the asymptotic form of e'^'^^^^ is not always true 
for the CME; note that there are chaotic behavior as well involved. If one considers a 
CME whose corresponding ODE is a 3-dimensional chaotic dynamics with a strange 
attractor, what will be the stationary distribution in the limit of V —>■ oo? This 
problem has been discussed in the past [19j. The general feeling is that is 
not smooth itself. So one does not have a "nice" 0(x)! For a very "rugged 0(x)", 
we believe that its Fenchel-Legendre transform c(A) might be a very powerful way 
to "find the key feature" of the The number of non-differentiable point is 

definitely much smaller than the number of peaks! 

Beyond deterministic dynamics. It is generally believed that when a system's 
size increases, the stochastic behavior at a mesoscopic level averaged out, and a 
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deterministic behavior emerges. However, our present analysis clearly show that 
the emerging deterministic behavior in the CME is a metastable system's dynamics. 
Beyond that time scale, another "macroscopic" stochastic behaviour exists! This 
multi-attractor stochastic system is a true emerging phenomenon that one can not 
naively expect from the deterministic dynamics (e.g., based on the relative area of the 
attractive basins) without detailed stochastic mechanistic modeling. The Maxwell 
construction is the consequence of the steady state on this "beyond-deterministic- 
infinite" time scale. 

There are three time scales in this mathematical hierarchy of cellular dynamics: 
A molecular signaling time scale (i.e., the rate constant for molecular interactions), 
a biochemical network time scale (i.e., the deterministic relaxation times to attrac- 
tors), and a cellular evolutionary time scale). We believe it is at the last level of 
stochastic dynamics that is most relevant to major cellular biological issues such as 
differentiation, apoptosis, and cancer immunoediting. 
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Figure 4: Schematics showing the mathematical hierarchy of cellular dynamics based 
on the chemical master equation (CME) approach, (a) stochastic dynamics based 
on the Gillespie algorithm; (b) deterministic dynamics tending to attractors; (c) 
probabilistic distributions for the two attractors; (d) stochastic dynamics among the 
attractors. (a), (b) and (c) represent stochastic molecular signaling, deterministic 
biochemical dynamics, and stochastic cellular evolution, respectively. 
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